/*
 * topological-hex, a program to compute combinatorial hexahedral meshes.
 *
 * Copyright (C) <2018> <Université catholique de Louvain (UCL), Belgique>
 *
 * List of the contributors to the development, description and complete
 * License: see LICENSE file.
 *
 * This program (topological-hex) is free software:
 * you can redistribute it and/or modify it under the terms
 * of the GNU General Public License as published by the Free
 * Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program (see COPYING file).  If not,
 * see <http://www.gnu.org/licenses/>.
 */

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <pthread.h>
#include <stdatomic.h>
#include <time.h>

#include "parallel_dfs.h"
#include "server.h"
#include "compatibility.h"

// #define MIN_VERTEX

typedef struct options {
  size_t num_subproblems;
  size_t chunk_size;

  size_t first_it;
  size_t max_it;

  hex_index    max_hex;
  vertex_index min_v_int;
  hex_index    min_removed_hex;
  vertex_index min_removed_v;

  const char *input_path;
  const char *output_prefix;

  const char *backup_dir;
  bool restore;
} options;

static void parse_options(options *opts, int argc, char **argv);

typedef struct cb_args {
  pthread_mutex_t *lock;
  quad_descriptor *descriptors;
  vertex_index num_vertices;

  uint32_t min_size;
  hex_index num_hexes;
  vertex_index *vertices;

  const hex_index *selected;
  hex_index num_selected;

  mesh *old_mesh;
  vertex_index *local_to_global;

  bool wrote_progress;
  double start_time;
} cb_args;

typedef struct on_listen_args {
  atomic_bool server_ready;
  pthread_cond_t cond;
  pthread_mutex_t lock;
} on_listen_args;

typedef struct server_args {
  cb_args *cb_args;
  on_listen_args *listen_args;
  server_backup_info *backup;
} server_args;

static int compare_ids(const void *ptr_a, const void *ptr_b) {
  hex_index a = *(const hex_index*)ptr_a;
  hex_index b = *(const hex_index*)ptr_b;

  if (a > b) return 1;
  if (a < b) return -1;
  return 0;
}

static double get_time(void) {
#ifdef SET_GLIBC_LINK_VERSIONS_HEADER
  return time(NULL);
#else
  struct timespec time;
  clock_gettime(CLOCK_MONOTONIC, &time);
  return time.tv_sec + 1e-9 * time.tv_nsec;
#endif
}

static bool array_contains(const hex_index *ids, hex_index size,
                           hex_index key);

static
void on_solution(
  const char *user_name,
  const vertex_index *vertices, hex_index num_hexes,
  void *data);

static void on_log(const char *message, void *data);
static void on_listen(void *data);
static void on_progress(uint32_t i, uint32_t n, void *data);

static
void *server_thread(void *data);

static void compute_adjacency(mesh *mesh);

static bool load_submesh(const mesh *mesh, hex_index *hex_ids,
                         hex_index *num_selected,
                         const char *backup_dir);

static void select_submesh(const mesh *mesh, hex_index *hex_ids,
                           hex_index *num_selected);

static void save_submesh(const char *backup_dir,
                         hex_index *ids, hex_index num_selected);

static void compute_boundary_vertices(
  const mesh *mesh,
  uint64_t *boundary_set, uint64_t *interior_set,
  const hex_index *hex_ids, hex_index num_selected);

static void build_vertex_map(const uint64_t *boundary_set,
                             const uint64_t *interior_set,
                             vertex_index *global_to_local,
                             vertex_index *local_to_global,
                             vertex_index *num_vertices);

static void build_quad_descriptors(
  const mesh *mesh,
  const hex_index *hex_ids, hex_index num_selected,
  const vertex_index *global_to_local, vertex_index num_vertices,
  quad_descriptor *local_descriptors,
  quad_descriptor *global_descriptors);

static void compute_boundary_quads(
  const mesh *mesh,
  const hex_index *hex_ids, hex_index num_selected,
  const vertex_index *global_to_local,
  vertex_index *quads, quad_index *num_quads);

static void show_cavity(mesh *mesh, hex_index *cavity, hex_index cavity_size,
                        uint64_t *interior_vertices);

static void setup_cavity_solver(solver *solver,
                                mesh *mesh, struct mesh *local_mesh,
                                hex_index *cavity, hex_index cavity_size,
                                uint64_t *boundary_vertices,
                                vertex_index *global_to_local,
                                vertex_index *local_to_global,
                                vertex_index min_removed_v);

static void remesh_cavity(mesh *mesh,
                          hex_index *old_cavity_hexes, hex_index old_cavity_size,
                          vertex_index *new_cavity, hex_index new_cavity_size);



int main(int argc, char **argv) {
  int status = 0;

  options opts;
  parse_options(&opts, argc, argv);

  mesh mesh;
  if (mesh_load(&mesh, opts.input_path) != SUCCESS) {
    status = 1;
    fprintf(stderr, "Failed to load mesh from file: %s\n",
            opts.input_path);
    goto fail_mesh_load;
  }

  pthread_mutex_t lock;
  if (pthread_mutex_init(&lock, NULL) != 0) {
    status = 1;
    fprintf(stderr, "Failed to initialize lock!\n");
    goto fail_init_lock;
  }

  compute_adjacency(&mesh);
//判断任意两个六面体的相容性
  for (hex_index i = 0; i < mesh.num_hexes; i++) {
    for (hex_index j = 0; j < i; j++) {
      if (!hex_compatibility(mesh.hexes + 8 * i, mesh.hexes + 8 * j))
        fprintf(stderr, "incompatibility: (%u, %u)\n", i, j);
    }
  }

  size_t iteration_id = opts.first_it;//迭代次数

  while (iteration_id < opts.max_it) {//开始迭代
    hex_index ids[64];
    hex_index num_selected;
    if (!opts.backup_dir ||
        !opts.restore ||
        !load_submesh(&mesh, ids, &num_selected, opts.backup_dir)) {
      select_submesh(&mesh, ids, &num_selected);
      opts.restore = false;
    }

    if (opts.backup_dir)
      save_submesh(opts.backup_dir, ids, num_selected);

    uint64_t boundary_vertices[2];
    uint64_t interior_vertices[2];
    compute_boundary_vertices(
      &mesh, boundary_vertices, interior_vertices,
      ids, num_selected);

    vertex_index local_to_global[128], global_to_local[128];
    vertex_index num_vertices;
    build_vertex_map(boundary_vertices, interior_vertices,
                     global_to_local, local_to_global,
                     &num_vertices);

    quad_descriptor local_descriptors[(uint32_t)num_vertices *
                                      (uint32_t)num_vertices];
    quad_descriptor global_descriptors[mesh.num_vertices * mesh.num_vertices];
    build_quad_descriptors(&mesh, ids, num_selected,
                           global_to_local, num_vertices,
                           local_descriptors, global_descriptors);

    struct mesh local_mesh;
    mesh_init(&local_mesh);

    local_mesh.num_vertices = num_vertices;
    posix_memalign((void**)&local_mesh.vertices, ALIGNMENT,
                   4 * num_vertices * sizeof(vertex_index));
    memset(local_mesh.vertices, 0, 4 * num_vertices * sizeof(vertex_index));

    posix_memalign((void**)&local_mesh.quads, ALIGNMENT,
                   4 * num_selected * 6 * sizeof(vertex_index));
    quad_index num_quads;
    compute_boundary_quads(&mesh, ids, num_selected, global_to_local,
                           local_mesh.quads, &num_quads);
    local_mesh.num_boundary_quads = local_mesh.num_quads = num_quads;

    if (num_selected <= opts.max_hex &&
        (popcount64(interior_vertices[0]) +
         popcount64(interior_vertices[1])) >= opts.min_v_int) {
      mesh_reserve(&local_mesh, num_selected - opts.min_removed_hex);

      show_cavity(&mesh, ids, num_selected, interior_vertices);
      vertex_index solution[8 * num_selected];

      cb_args args;
      args.lock = &lock;
      args.descriptors = local_descriptors;
      args.num_vertices = num_vertices;
#ifdef MIN_VERTEX
      args.min_size = num_vertices;
#else
      args.min_size = num_selected;
#endif
      args.num_hexes = num_selected;
      args.vertices = solution;
      args.local_to_global = local_to_global;
      args.old_mesh = &mesh;
      args.num_selected = num_selected;
      args.selected = ids;
      args.start_time = get_time();
      args.wrote_progress = false;

      uint32_t start_size = args.min_size;

      solver solver;
      setup_cavity_solver(&solver, &mesh, &local_mesh,
                          ids, num_selected, boundary_vertices,
                          global_to_local, local_to_global,
                          opts.min_removed_v);

      on_listen_args listen_args;
      atomic_init(&listen_args.server_ready, false);
      pthread_mutex_init(&listen_args.lock, NULL);
      pthread_cond_init(&listen_args.cond, NULL);

      server_backup_info info = {opts.backup_dir, opts.restore};
      server_args server_args = {&args, &listen_args, NULL};
      if (opts.backup_dir) {
        server_args.backup = &info;
      }

      pthread_t thread;
      pthread_create(&thread, NULL, server_thread, &server_args);

      pthread_mutex_lock(&listen_args.lock);
      while (!atomic_load_explicit(&listen_args.server_ready,
                                   memory_order_acquire)) {
        pthread_cond_wait(&listen_args.cond, &listen_args.lock);
      }
      pthread_mutex_unlock(&listen_args.lock);

      if (opts.restore) {
        solver_run_distributed_client("localhost", on_log, &args,
                                      opts.chunk_size);
      }
      else {
        solver_run_distributed(&solver, "localhost",
                               on_log, &args,
                               opts.chunk_size,
                               opts.num_subproblems);
      }

      solver_release(&solver);

      opts.restore = false;

      pthread_join(thread, NULL);

      pthread_cond_destroy(&listen_args.cond);
      pthread_mutex_destroy(&listen_args.lock);

      printf("\n");

      if (args.min_size < start_size) {
        for (uint32_t i = 0; i < 8 * args.num_hexes; i++)
          solution[i] = local_to_global[solution[i]];

        remesh_cavity(&mesh, ids, num_selected, solution, args.num_hexes);

        printf("Reduced mesh size to %u\n", mesh.num_hexes);

        char filename[256];
        snprintf(filename, sizeof(filename), "%s-%zu.mesh",
                 opts.output_prefix, iteration_id++);

        save_mesh(&mesh, filename);
      }
    }

    mesh_release(&local_mesh);
  }

  pthread_mutex_destroy(&lock);
fail_init_lock:
  mesh_release(&mesh);
fail_mesh_load:
  return status;
}

static void print_usage(const char *program_name);

static size_t process_integer(int *argc, char **argv, int i,
                              size_t min, size_t max);

static void missing_argument(const char *program_name, const char *option);
static void bad_integer(const char *program_name,
                        const char *option, const char *arg);

static void parse_options(options *opts, int argc, char **argv) {
  opts->num_subproblems = 1024 * 1024;
  opts->chunk_size = 64;

  opts->first_it = 0;
  opts->max_it = 256 * 1024;

  opts->max_hex = 10;
  opts->min_v_int = 3;
  opts->min_removed_hex = 2;
  opts->min_removed_v = 1;

  opts->output_prefix = "result";
  opts->backup_dir = NULL;
  opts->restore = false;

  int i = 1;
  while (i < argc) {
    if (strcmp(argv[i], "-o") == 0 || strcmp(argv[i], "--output") == 0) {
      if (i + 1 == argc) missing_argument(argv[0], argv[i]);

      opts->output_prefix = argv[i + 1];
      memmove(argv + i, argv + i + 2, (argc - i - 2) * sizeof(*argv));

      argc -= 2;
    }
    else if (strcmp(argv[i], "--chunk-size") == 0)
      opts->chunk_size = process_integer(&argc, argv, i, 1, 1024 * 1024);
    else if (strcmp(argv[i], "--num-subproblems") == 0)
      opts->num_subproblems = process_integer(&argc, argv, i, 1, SIZE_MAX);
    else if (strcmp(argv[i], "--first-it") == 0)
      opts->first_it = process_integer(&argc, argv, i, 0, SIZE_MAX);
    else if (strcmp(argv[i], "--max-it") == 0)
      opts->max_it = process_integer(&argc, argv, i, 0, SIZE_MAX);
    else if (strcmp(argv[i], "--max-hex") == 0)
      opts->max_hex = process_integer(&argc, argv, i, 0, 256);
    else if (strcmp(argv[i], "--min-v-int") == 0)
      opts->min_v_int = process_integer(&argc, argv, i, 0, 256);
    else if (strcmp(argv[i], "--min-v-int") == 0)
      opts->min_v_int = process_integer(&argc, argv, i, 0, 256);
    else if (strcmp(argv[i], "--min-removed-hex") == 0)
      opts->min_removed_hex = process_integer(&argc, argv, i, 0, 256);
    else if (strcmp(argv[i], "--min-removed-v") == 0)
      opts->min_removed_v = process_integer(&argc, argv, i, 0, 256);
    else if (strcmp(argv[i], "--backup") == 0) {
      if (i + 1 == argc) missing_argument(argv[0], argv[i]);

      opts->backup_dir = argv[i + 1];
      memmove(argv + i, argv + i + 2, (argc - i - 2) * sizeof(*argv));

      argc -= 2;
    }
    else if (strcmp(argv[i], "--restore") == 0) {
      opts->restore = true;
      memmove(argv + i, argv + i + 1, (argc - i - 1) * sizeof(*argv));
      argc--;
    }
    else if (strcmp(argv[i], "--seed") == 0)
      srand(process_integer(&argc, argv, i, 0, UINT32_MAX));
    else
      i++;
  }

  if (argc == 2)
    opts->input_path = argv[1];
  else if (argc > 2) {
    fprintf(stderr, "%s: extranumerous arguments after input\n", argv[0]);
    print_usage(argv[0]);
    exit(1);
  }
  else {
    fprintf(stderr, "%s: missing input file\n", argv[0]);
    print_usage(argv[0]);
    exit(1);
  }
}

static void print_usage(const char *program_name) {
  printf("Usage: %s [OPTIONS] INPUT_MESH\n", program_name);
  printf("\n");
  printf("  -o, --output PREFIX\n");
  printf("    Specifies a prefix for output files (default: result)\n");
  printf("  --chunk-size N\n");
  printf("    Specifies the minimum unit of work done by the local solver.\n");
  printf("  --num-suproblems N\n");
  printf("    Specifies how many subproblems should be created for parallel\n");
  printf("    exploration.\n");
  printf("  --first-it N\n");
  printf("    Specifies the index of the first iteration. Used to format output "
         "files.\n");
  printf("  --max-it N\n");
  printf("    Maximum number of iteration to perform.\n");
  printf("  --max-hex H\n");
  printf("    Maximum number hexahedra in a cavity.\n");
  printf("  --min-v-int V\n");
  printf("    Minimum number of interior vertices in a cavity.\n");
  printf("  --min-removed-hex H\n");
  printf("    Minimum number of hexahedra removed by an iteration.\n");
  printf("  --min-removed-v V\n");
  printf("    Minimum number of vertices removed by an iteration.\n");
  printf("  --backup DIRECTORY\n");
  printf("    Directory used to store backups if the program is interrupted during");
  printf("    its execution.\n");
  printf("  --backup DIRECTORY\n");
  printf("    Directory used to store backup.\n");
  printf("  --restore\n");
  printf("    Restores a backup from the backup directory.\n");
  printf("  --seed SEED\n");
  printf("    Integer seed for the random number generator.\n");
}

static size_t process_integer(int *argc, char **argv, int i,
                              size_t min, size_t max) {
  if (i + 1 == *argc) missing_argument(argv[0], argv[i]);

  char *end = NULL;
  size_t result = strtoull(argv[i + 1], &end, 10);
  if (end == NULL || *end != '\0')
    bad_integer(argv[0], argv[i], argv[i + 1]);

  if (result < min || result >= max) {
    fprintf(stderr, "%s: out of range argument for option `%s': "
            "%zu (min: %zu, max: %zu)\n", argv[0], argv[i], result,
            min, max);
  }

  memmove(argv + i, argv + i + 2, (*argc - i - 2) * sizeof(*argv));
  *argc -= 2;

  return result;
}

static void missing_argument(const char *program_name, const char *option) {
  fprintf(stderr, "%s: missing argument for option `%s'\n", program_name, option);
  print_usage(program_name);
  exit(1);
}

static void bad_integer(const char *program_name,
                        const char *option, const char *arg) {
  fprintf(stderr, "%s: bad format for option `%s': `%s'\n",
          program_name, option, arg);
  print_usage(program_name);
  exit(1);
}

static void on_log(const char *message, void *data) {
  /* cb_args *args = data; */
  /* pthread_mutex_lock(args->lock); */
  /* printf("%s\n", message); */
  /* pthread_mutex_unlock(args->lock); */
}

static void on_listen(void *data) {
  on_listen_args *args = data;
  atomic_store_explicit(&args->server_ready, true, memory_order_release);
  pthread_cond_broadcast(&args->cond);
}

#define BAR_WIDTH 50

static void on_progress(uint32_t i, uint32_t n, void *data) {
  cb_args *args = data;

  FILE *out = stdout;

  pthread_mutex_lock(args->lock);
  if (!args->wrote_progress)
    fprintf(out, "\n");

  char bar[BAR_WIDTH + 1];
  uint64_t ticks = (BAR_WIDTH * (uint64_t)i) / n;
  memset(bar, '=', ticks);
  if (ticks != 0 && i != n) {
    bar[ticks - 1] = '>';
    memset(bar + ticks, ' ', BAR_WIDTH - ticks);
  }
  else if (ticks == 0)
    memset(bar + ticks, ' ', BAR_WIDTH - ticks);
  bar[BAR_WIDTH] = 0;

  fprintf(out, "\r\033[K[%s] ", bar);

  if (i == n)
    fprintf(out, "DONE!");
  else if (i != 0) {
    double elapsed = get_time() - args->start_time;
    uint64_t eta = (n - i) * (elapsed / i);

    uint64_t days = eta / (24 * 3600);
    eta -= 24 * 3600 * days;
    uint64_t hours = eta / 3600;
    eta -= 3600 * hours;
    uint64_t minutes = eta / 60;
    uint64_t seconds = eta - 60 * minutes;

    if (days != 0)
      fprintf(out, "ETA: %zud%02zu:%02zu:%02zu", days, hours, minutes, seconds);
    else if (hours != 0)
      fprintf(out, "ETA: %02zu:%02zu:%02zu", hours, minutes, seconds);
    else if (minutes != 0)
      fprintf(out, "ETA: %02zu:%02zu", minutes, seconds);
    else
      fprintf(out, "ETA: %zus", seconds);
  }

  args->wrote_progress = true;

  fflush(out);
  pthread_mutex_unlock(args->lock);
}

static
void *server_thread(void *data) {
  server_args *args = data;
  if (server_run(on_solution, args->cb_args,
                 on_log, args->cb_args,
                 on_listen, args->listen_args,
                 on_progress, args->cb_args,
                 args->backup) != SUCCESS) {
    fprintf(stderr, "Failed to start server!\n");
    on_listen(args->listen_args);
  }
  return NULL;
}

static
void on_solution(const char *user_name,
                 const vertex_index *vertices, hex_index num_hexes,
                 void *data) {
  cb_args *args = data;

  pthread_mutex_lock(args->lock);
  if (args->wrote_progress)
    printf("\n");
  args->wrote_progress = false;

  printf("SOLUTION from %s!\n", user_name);

  bool compatible = true;

  for (hex_index i = 0; i < num_hexes; i++) {
    vertex_index global_hex[8];
    for (hex_vertex_index k = 0; k < 8; k++)
      global_hex[k] = args->local_to_global[vertices[8*i + k]];

    for (hex_index j = 0; j < i; j++) {
      if (!hex_compatibility(vertices + 8*i, vertices + 8*j))
        printf("Incompatibility: %u, %u\n", i, j);
    }

    for (hex_index j = 0; compatible && j < args->old_mesh->num_hexes; j++) {
      if (array_contains(args->selected, args->num_selected, j))
        continue;

      if (!hex_compatibility(global_hex, args->old_mesh->hexes + 8*j)) {
        printf("  Hex #%u is incompatible with external hex #%u!\n",
               i, j);
        printf("    A: {");
        for (hex_vertex_index k = 0; k < 8; k++) {
          printf("%u", global_hex[k]);
          if (k + 1 != 8) printf(", ");
        }
        printf("}\n");
        printf("    B: {");
        for (hex_vertex_index k = 0; k < 8; k++) {
          printf("%u", args->old_mesh->hexes[8*j + k]);
          if (k + 1 != 8) printf(", ");
        }
        printf("}\n");
        compatible = false;
      }
    }

    for (hex_facet_index j = 0; j < 6; j++) {
      const hex_vertex_index *indices = HEX_FACET_CONTENTS[j];

      vertex_index quad[4];
      for (quad_vertex_index k = 0; k < 4; k++)
        quad[k] = vertices[8*i + indices[k]];

      quad_descriptor desc = args->descriptors[
        diagonal_index(args->num_vertices, quad[0], quad[2])];
      hex_index adjacent = quad_descriptor_hex_index(desc);
      if (desc != EMPTY_QUAD_DESCRIPTOR && adjacent != NO_HEXAHEDRON) {
        for (hex_facet_index k = j + 1; k < 6; k++) {
          const hex_vertex_index *other_indices = HEX_FACET_CONTENTS[k];

          vertex_index other_quad[4];
          for (quad_vertex_index l = 0; l < 4; l++)
            other_quad[l] = vertices[8*i + other_indices[l]];

          quad_descriptor other_desc = args->descriptors[
            diagonal_index(args->num_vertices, other_quad[0], other_quad[2])];
          hex_index other_adjacent = quad_descriptor_hex_index(other_desc);
          if (other_desc != EMPTY_QUAD_DESCRIPTOR && adjacent == other_adjacent) {
            printf("  Hex #%u shares two faces with external hex #%u!\n",
                   i, adjacent);
            compatible = false;
          }
        }
      }
    }

    printf("  %u) {", i);
    for (hex_vertex_index j = 0; j < 8; j++) {
      printf("%u", vertices[8*i+j]);
      if (j + 1 != 8) printf(", ");
    }
    printf("}\n");
  }

  uint32_t size = num_hexes;
#ifdef MIN_VERTEX
  size = 0;
  for (uint32_t i = 0; i < 8 * num_hexes; i++) {
    if (size < vertices[i] + 1u)
      size = vertices[i] + 1u;
  }
#endif

  printf("new size: %u (vs. %u)\n", size, args->min_size);

  if (compatible &&
      (size < args->min_size ||
       (size == args->min_size && num_hexes < args->num_hexes))) {
    args->min_size = size;
    args->num_hexes = num_hexes;
    memcpy(args->vertices, vertices,
           8 * num_hexes * sizeof(*args->vertices));
  }

  pthread_mutex_unlock(args->lock);
}
//计算六面体周围的面和六面体
static void compute_adjacency(mesh *mesh) {
  quad_descriptor descriptors[mesh->num_vertices * mesh->num_vertices];
  for (uint32_t i = 0; i < mesh->num_vertices * mesh->num_vertices; i++)
    descriptors[i] = EMPTY_QUAD_DESCRIPTOR;

  for (hex_index hex_id = 0; hex_id < mesh->num_hexes; hex_id++) {
    vertex_index *hex = mesh->hexes + 8 * hex_id;

    for (hex_facet_index i = 0; i < 6; i++) {
      const hex_vertex_index *indices = HEX_FACET_CONTENTS[i];

      vertex_index quad[4];
      for (quad_vertex_index j = 0; j < 4; j++)
        quad[j] = hex[indices[j]];

      quad_descriptor old_descriptor =
        descriptors[diagonal_index(mesh->num_vertices, quad[0], quad[2])];
      if (old_descriptor != EMPTY_QUAD_DESCRIPTOR) {
        mesh->hex_adjacency[6*hex_id + i] =
          quad_descriptor_hex_index(old_descriptor);
        mesh->hex_facets[6*hex_id + i] =
          quad_descriptor_quad_index(old_descriptor);

        mesh->hex_adjacency[6*quad_descriptor_hex_index(old_descriptor) +
                            quad_descriptor_quad_index(old_descriptor)] =
          hex_id;
        mesh->hex_facets[6*quad_descriptor_hex_index(old_descriptor) +
                         quad_descriptor_quad_index(old_descriptor)] =
          i;
      }
      else {
        descriptors[diagonal_index(mesh->num_vertices, quad[0], quad[2])] =
          make_quad_descriptor(i, hex_id, quad[1], quad[3]);

        descriptors[diagonal_index(mesh->num_vertices, quad[1], quad[3])] =
          make_quad_descriptor(i, hex_id, quad[0], quad[2]);

        mesh->hex_adjacency[6*hex_id + i] = NO_HEXAHEDRON;
        mesh->hex_facets[6*hex_id + i] = 0;
      }
    }
  }
}

static bool load_submesh(const mesh *mesh, hex_index *hex_ids,
                         hex_index *num_selected,
                         const char *backup_dir) {
  char filename[1024];
  snprintf(filename, sizeof(filename), "%s/cavity.data", backup_dir);//读入文件在当前目录下，名为cavity.data

  FILE *file = fopen(filename, "r");
  if (!file)
    return false;

  bool success = true;

  if (fread(num_selected, sizeof(*num_selected), 1, file) == 1) {//读取num_selected和hex_ids
    success = (fread(hex_ids, sizeof(*hex_ids) * (*num_selected), 1,
                     file) == 1);
  }
  else
    success = false;

  fclose(file);
  return success;
}

static void select_submesh(const mesh *mesh,
                           hex_index *hex_ids, hex_index *num_selected) {
  uint64_t selected_hex[2] = {0};

  hex_index start_id = rand() % mesh->num_hexes;//start_id随机选择为0到num_hex-1的某一个

  *num_selected = 0;//选择六面体数量初始化为0

  hex_index queue[mesh->num_hexes * 6];//队列最大为mesh中六面体数量的6倍
  size_t queue_size = 0, queue_offset = 0;

  queue[queue_size++] = start_id;//队列第一个元素是start_id

  hex_index max = 8 + (rand() % 20);//max为8到27中的随机数

  while (*num_selected < max && queue_offset != queue_size) {//初始化时queue_size是1，queue_offset是0
    hex_index id = queue[queue_offset++];//id是队尾元素，首次循环为start_id
    if (queue_size >= sizeof(queue)/sizeof(*queue) - 6) {//队列元素数量和允许最大数量相差6时
      memmove(queue, queue + queue_offset,
              sizeof(queue) - sizeof(*queue) * queue_offset);//把offset后的部分copy到队首，覆盖队首内容
      queue_size -= queue_offset;//队列清空
      queue_offset = 0;
    }
    //只考虑id为0-63的情况，id/64=0,id%64=id;
    if ((selected_hex[id / 64] & (1ull << (id % 64))) == 0) {//首次循环由于selected_hex都是0，恒成立,第一次循环s_h[0]=2^id,运算结果为2^id，不成立
      hex_ids[(*num_selected)++] = id;//把id记录的六面体放入选择数组中
      //若id为0到63，s_h[0] |= 2^id,首次循环时，s_h[0]赋值2^id
      //第一次循环，s_h[0]=2^id,s_h[0]|=2^id,结果不变
      selected_hex[id / 64] |= (1ull << (id % 64));
    }
    //从队尾六面体随机选择0-6个方向的六面体入队
    for (hex_facet_index i = 0; i < 6; i++) {
      if (mesh->hex_adjacency[6*id + i] != NO_HEXAHEDRON &&
          queue_size != sizeof(queue)/sizeof(*queue)) {
        if (rand() % 2) queue[queue_size++] = mesh->hex_adjacency[6*id + i];
      }
    }
  }

  qsort(hex_ids, *num_selected, sizeof(*hex_ids), compare_ids);
}

static void save_submesh(const char *backup_dir,
                         hex_index *ids, hex_index num_selected) {
  char filename[1024];
  snprintf(filename, sizeof(filename), "%s/cavity.data", backup_dir);

  FILE *file = fopen(filename, "w");
  if (!file)
    return;

  fwrite(&num_selected, sizeof(num_selected), 1, file);
  fwrite(ids, sizeof(*ids) * num_selected, 1, file);
  fclose(file);
}

static bool array_contains(const hex_index *ids, hex_index size,
                           hex_index key) {
  for (hex_index i = 0; i < size; i++)
    if (ids[i] == key) return true;
  return false;
}

static void compute_boundary_vertices(
  const mesh *mesh,
  uint64_t *boundary_set, uint64_t *interior_set,
  const hex_index *hex_ids, hex_index num_selected) {
  boundary_set[0] = boundary_set[1] = 0;
  interior_set[0] = interior_set[1] = 0;

  for (hex_index i = 0; i < num_selected; i++) {
    hex_index hex_id = hex_ids[i];

    for (hex_facet_index j = 0; j < 6; j++) {
      if (!array_contains(hex_ids, num_selected,
                          mesh->hex_adjacency[6*hex_id + j])) {
        for (quad_vertex_index k = 0; k < 4; k++) {
          vertex_index v = mesh->hexes[8*hex_id + HEX_FACET_CONTENTS[j][k]];
          boundary_set[v / 64] |= (1ull << (v % 64));
        }
      }
    }
  }

  for (hex_index i = 0; i < num_selected; i++) {
    hex_index hex_id = hex_ids[i];

    for (hex_vertex_index j = 0; j < 8; j++) {
      vertex_index v = mesh->hexes[8*hex_id + j];
      if ((boundary_set[v / 64] & (1ull << (v % 64))) == 0)
        interior_set[v / 64] |= (1ull << (v % 64));
    }
  }
}

static void build_vertex_map(const uint64_t *boundary_set,
                             const uint64_t *interior_set,
                             vertex_index *global_to_local,
                             vertex_index *local_to_global,
                             vertex_index *num_vertices) {
  vertex_index local_id = 0;

  for (vertex_index global_id = 0; global_id < 128; global_id++) {
    if (boundary_set[global_id / 64] & (1ul << (global_id % 64))) {
      global_to_local[global_id] = local_id;
      local_to_global[local_id]  = global_id;

      local_id++;
    }
  }

  for (vertex_index global_id = 0; global_id < 128; global_id++) {
    if (interior_set[global_id / 64] & (1ul << (global_id % 64))) {
      global_to_local[global_id] = local_id;
      local_to_global[local_id]  = global_id;

      local_id++;
    }
  }

  *num_vertices = local_id;
}

static void build_quad_descriptors(
  const mesh *mesh,
  const hex_index *hex_ids, hex_index num_selected,
  const vertex_index *global_to_local, vertex_index num_vertices,
  quad_descriptor *local_descriptors,
  quad_descriptor *global_descriptors) {
  for (uint32_t i = 0; i < (uint32_t)num_vertices * (uint32_t)num_vertices; i++)
    local_descriptors[i] = EMPTY_QUAD_DESCRIPTOR;

  for (uint32_t i = 0; i < mesh->num_vertices * mesh->num_vertices; i++)
    global_descriptors[i] = EMPTY_QUAD_DESCRIPTOR;

  for (hex_index i = 0; i < num_selected; i++) {
    hex_index hex_id = hex_ids[i];
    vertex_index *hex = mesh->hexes + 8 * hex_id;

    for (hex_facet_index j = 0; j < 6; j++) {
      if (!array_contains(hex_ids, num_selected,
                          mesh->hex_adjacency[6*hex_id + j])) {
        const hex_vertex_index *indices = HEX_FACET_CONTENTS[j];

        vertex_index quad[4];
        for (quad_vertex_index k = 0; k < 4; k++)
          quad[k] = hex[indices[k]];

        global_descriptors[
          diagonal_index(mesh->num_vertices, quad[0], quad[2])] =
          make_quad_descriptor(0, mesh->hex_adjacency[6*hex_id + j],
                               quad[1], quad[3]);

        global_descriptors[
          diagonal_index(mesh->num_vertices, quad[1], quad[3])] =
          make_quad_descriptor(0, mesh->hex_adjacency[6*hex_id + j],
                               quad[0], quad[2]);

        local_descriptors[
          diagonal_index(num_vertices,
                         global_to_local[quad[1]],
                         global_to_local[quad[3]])] =
          make_quad_descriptor(
            0, mesh->hex_adjacency[6*hex_id + j],
            global_to_local[quad[0]], global_to_local[quad[2]]);

        local_descriptors[
          diagonal_index(num_vertices,
                         global_to_local[quad[0]],
                         global_to_local[quad[2]])] =
          make_quad_descriptor(
            0, mesh->hex_adjacency[6*hex_id + j],
            global_to_local[quad[1]], global_to_local[quad[3]]);
      }
    }
  }
}

static void compute_boundary_quads(
  const mesh *mesh,
  const hex_index *hex_ids, hex_index num_selected,
  const vertex_index *global_to_local,
  vertex_index *quads, quad_index *num_quads) {
  *num_quads = 0;
  for (hex_index i = 0; i < num_selected; i++) {
    hex_index hex_id = hex_ids[i];
    vertex_index *hexes = mesh->hexes + 8 * hex_id;

    for (hex_facet_index j = 0; j < 6; j++) {
      if (!array_contains(hex_ids, num_selected,
                          mesh->hex_adjacency[6*hex_id + j])) {
        quad_vertex_index perm[] = {0, 3, 2, 1};
        for (quad_vertex_index k = 0; k < 4; k++) {
          quads[4 * (*num_quads) + k] =
            global_to_local[hexes[HEX_ORIENTED_FACET_CONTENTS[j][perm[k]]]];
        }

        (*num_quads)++;
      }
    }
  }
}

static void show_cavity(mesh *mesh, hex_index *cavity, hex_index cavity_size,
                        uint64_t *interior_vertices) {
  printf("selected: ");
  for (hex_index i = 0; i < cavity_size; i++)
    printf("%u ", cavity[i]);
  printf("\n");

  for (hex_index i = 0; i < cavity_size; i++) {
    for (hex_vertex_index j = 0; j < 8; j++) {
      printf("%u ", 1 + mesh->hexes[8*cavity[i] + j]);
    }
    printf("0\n");
  }

  printf("interior vertices: ");
  for (size_t i = 0; i < 128; i++) {
    if (interior_vertices[i / 64] & (1ul << (i % 64))) {
      printf("%zu ", i);
    }
  }
  printf("\n");

  printf("num selected: %u\n", cavity_size);
  printf("num quads: %u\n", mesh->num_quads);
}

static void setup_cavity_solver(solver *solver,
                                mesh *mesh, struct mesh *local_mesh,
                                hex_index *cavity, hex_index cavity_size,
                                uint64_t *boundary_vertices,
                                vertex_index *global_to_local,
                                vertex_index *local_to_global,
                                vertex_index min_removed_v) {
  solver_init(solver, local_mesh, local_mesh->num_vertices - min_removed_v,
              SOLVER_BIPARTITE_GRAPH);

  for (hex_index i = 0; i < mesh->num_hexes; i++) {
    if (array_contains(cavity, cavity_size, i)) continue;

    vertex_index *hex = mesh->hexes + 8 * i;

    for (hex_facet_index j = 0; j < 6; j++) {
      vertex_index global_face[4];

      for (quad_vertex_index k = 0; k < 4; k++)
        global_face[k] = hex[HEX_FACET_CONTENTS[j][k]];

      for (quad_vertex_index k = 0; k < 4; k++) {
        if ((boundary_vertices[global_face[k] / 64] &
             (1ul << (global_face[k] % 64))) == 0)
          continue;

        vertex_index v = global_to_local[global_face[k]];

        if ((boundary_vertices[global_face[(k+1)%4] / 64] &
             (1ul << (global_face[(k+1)%4] % 64)))) {
          vertex_index v_next = global_to_local[global_face[(k+1) % 4]];

          vertex_set_insert(solver->known_neighbors + v, v_next);
          vertex_set_insert(solver->known_neighbors + v_next, v);
        }

        if ((boundary_vertices[global_face[(k+2)%4] / 64] &
             (1ul << (global_face[(k+2)%4] % 64)))) {
          vertex_index v_across = global_to_local[global_face[(k+2) % 4]];

          vertex_set_remove(solver->allowed_neighbors + v, v_across);
          vertex_set_remove(solver->allowed_neighbors + v_across, v);

          vertex_set_insert(solver->known_quad_diagonals + v, v_across);
          vertex_set_insert(solver->known_quad_diagonals + v_across, v);

          if ((boundary_vertices[global_face[(k+1)%4] / 64] &
               (1ul << (global_face[(k+1)%4] % 64))) == 0 ||
              (boundary_vertices[global_face[(k+3)%4] / 64] &
               (1ul << (global_face[(k+3)%4] % 64))) == 0) {
            solver->quad_descriptors[
              diagonal_index(solver->num_vertices, v, v_across)] =
              make_quad_descriptor(0, NO_HEXAHEDRON,
                                   VERTEX_SET_MAX, VERTEX_SET_MAX);
          }
        }
      }
    }

    for (quad_vertex_index k = 0; k < 4; k++) {
      if ((boundary_vertices[hex[k] / 64] &
           (1ul << (hex[k] % 64))) == 0)
        continue;

      if ((boundary_vertices[hex[4+((k+2)%4)] / 64] &
           (1ul << (hex[4+((k+2)%4)] % 64))) == 0)
        continue;

      vertex_index v = global_to_local[hex[k]];
      vertex_index v_across = global_to_local[hex[4 + ((k+2) % 4)]];

      vertex_set_remove(solver->allowed_neighbors + v, v_across);
      vertex_set_remove(solver->allowed_neighbors + v_across, v);

      vertex_set_insert(solver->known_diagonals + v, v_across);
      vertex_set_insert(solver->known_diagonals + v_across, v);
    }
  }
}

static void remesh_cavity(mesh *mesh,
                          hex_index *old_cavity_hexes, hex_index old_cavity_size,
                          vertex_index *new_cavity, hex_index new_cavity_size) {
  for (hex_index i = 0; i < new_cavity_size; i++) {
    memcpy(mesh->hexes + 8 * old_cavity_hexes[i],
           new_cavity + 8 * i,
           8 * sizeof(*new_cavity));
  }

  for (hex_index i = 0; i < old_cavity_size - new_cavity_size; i++) {
    hex_index to_remove = old_cavity_hexes[old_cavity_size - i - 1];
    memmove(mesh->hexes + 8 * to_remove,
            mesh->hexes + 8 * (to_remove + 1),
            8 * sizeof(*mesh->hexes) * (mesh->num_hexes - to_remove - 1));
    mesh->num_hexes--;
  }

  compute_adjacency(mesh);
}
